namespace Gsl
{
    using System;
    using System.Runtime.InteropServices;

    public class LeastSquaresFit
    {
        [DllImport("libgsl")]
        private static extern int gsl_fit_linear(double[] x, int xstride, double[] y, int ystride, int n, out double c0, out double c1, out double cov00, out double cov01, out double cov11, out double sumsq);
        // warning:  this function needs a function in libgslcblas, which
        // is not pulled in because of bad linking (I guess).
        // So if you need this function, do a
        // "LD_PRELOAD=/usr/lib/libgslcblas.so"
        // before running your program
        // or ... (new) ... just add the "sqrt" function in gslcblas and run it somewhere
        [DllImport("libgsl")]
        private static extern int gsl_multifit_linear(Matrix X, Vector y, Vector c, Matrix cov, out double chisq, IntPtr work);
        [DllImport("libgsl")]
        private static extern IntPtr gsl_multifit_linear_alloc(int n, int p);
        [DllImport("libgsl")]
        private static extern void gsl_multifit_linear_free(IntPtr w);
        [DllImport("libgslcblas")]
        private static extern double sqrt(double x); 
        
        public static FitInfo FitLinear(double[] x, double[] y)
        {
            FitInfo info = new FitInfo(2);
            
            double c0, c1;
            double cov00, cov01, cov11;            
            Gsl.Error.Code ret = (Gsl.Error.Code)gsl_fit_linear(x, 1, y, 1, x.Length, out c0, out c1, out cov00, out cov01, out cov11, out info.SumSquares);
            
            info.C[0] = c0;
            info.C[1] = c1;
            info.Covariance[0, 0] = cov00;
            info.Covariance[0, 1] = cov01;
            info.Covariance[1, 1] = cov11;
            
            if (ret != Gsl.Error.Code.GSL_SUCCESS)
                throw new GslException(String.Format("GSL error in linear fitting: {0}", ret), ret);
                
            return info;
        }
        
        public static FitInfo FitMultiparameter(Matrix x, Vector y)
        {
            FitInfo info = new FitInfo(x.NumColumns);
            
            sqrt(1.0);
            IntPtr w = gsl_multifit_linear_alloc(x.NumRows, x.NumColumns);
            
            Gsl.Error.Code ret = (Gsl.Error.Code)gsl_multifit_linear(x, y, info.C, info.Covariance, out info.SumSquares, w);
            
            gsl_multifit_linear_free(w);
            
            if (ret != Gsl.Error.Code.GSL_SUCCESS)
                throw new GslException(String.Format("GSL error in linear fitting: {0}", ret), ret);
            
            return info;
        }
        
        public static FitInfo FitPolynomial(int index, double[] x, double[] y)
        {
            Matrix xx = new Matrix(x.Length, index);
            Vector yy = new Vector(x.Length);
            for (int i = 0; i < x.Length; i++)
            {
                for (int j = 0; j < index; j++)
                    xx[i, j] = Math.Pow(x[i], (double)j);
                yy[i] = y[i];
            }
            
            return FitMultiparameter(xx, yy);
        }
    }
}
